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Abstract 

We investigate the emergence of long-range correlations in granular shear flow. By increasing 
the density of a simulated granular flow we observe a spontaneous transition from a dilute regime, 
where interactions are dominated by binary collisions, to a dense regime characterized by large 
force networks and collective motions. With increasing density, interacting grains tend to form 
networks of simultaneous contacts due to the dissipative nature of collisions. We quantify the size 
of these networks by measuring correlations between grain forces and find that there are dramatic 
changes in the statistics of contact forces as the size of the networks increases. 
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I. INTRODUCTION 



Granular materials exhibit a wide range of fascinating behaviors [1], but predictive the- 
ories linking the microscopic grain interactions to macroscopic properties remains a topic 
of much debate. Of particular interest are constitutive relations for granular flow, which 
are important for engineering and geophysical applications 

aa 

and challenge the tenets 

of conventional statistical physics [4J. Because granular materials are athermal their dy- 
namics always occur far from equilibrium and a proper formulation of constitutive relations 
relies on the construction of non-equilibrium statistical theories that must be sensitive to 
the interactions between grains. 

Interactions in realistic granular materials arise due to grain elasticity and friction, but are 
complicated by various other mechanisms including humidity la la |7|, |8|, |9j, grain shapes 
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12j, and fracture processes occurring within the material 13j. However, a great deal of 
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theoretical and computational progress has been made using the simple approximation that 

n 

grains are spherical and perfectly dry [141] . In this case a purely repulsive force arises when 
two grains come into contact due to the deformation of grains and friction between grains. 

Understanding the nature of grain forces and dynamics, even in this relatively simple 
case, has proven difficult. At very low densities it can safely be assumed that only binary 
interactions occur and constitutive relations can be determined by statistically tracking the 
repulsive force created in each interaction. This is the basis of kinetic theory, which has been 



successfully applied to granular flows [15. However, for very large densities, it is observed 



that multi-grain contacts always occur and contact forces are transmitted through 



21 



22 



23|]. For 



"force chain networks" formed by the topology of the contact network 
these high densities the forces between contacting grains still arise from grain deformation 
and friction, but the extent of the interactions is not localized and depends on properties of 
the force chain networks. 



The presence of force chain networks calls into question theories that assume loca 



interactions and has inspired new models based on properties of the force chains 24], 



3, 
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ized 
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301 ] . However, although force networks can be visualized, it has proven difficult to 



31 



32 



33 



341 ] . This has led to the 



measure quantitative correlations between contact forces 
speculation that force chain networks are simply a perceived correlation, until recently when 
long-range correlations were measured between the averaged contact forces in a quasi-static 
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experimental shear flow at high density 



23]. 



This discovery raises important questions about the proper assumptions to make when 
constructing theories of granular materials. For very dilute systems only binary collisions 
occur and force chain networks do not play a role; for very dense systems force networks are 
the dominating microscopic interaction. In order to understand the origin of macroscopic 
properties in granular flows, it is necessary to pinpoint exactly how and when correlations 
appear. 

In this paper we measure spatial correlations of the total force on grains undergoing shear 
deformation. This measurement defines a characteristic length-scale £ quantifying the size 
of force networks arising from clusters of simultaneously contacting grains. We find that £ 
grows with packing fraction, diverges at a finite packing fraction, and has measurable effects 
on the contact forces between grains. The correlation measurement also provides a natural 
boundary between dilute flows where only binary collisions occur and dense flows where 
force networks spontaneously emerge. 

We begin in Section [Til by briefly outlining the phenomenology of granular shear flow 
to define the exact regime we will be studying. In Section II III we discuss the numerical 
algorithm used to simulate shear flow. In Section HVl we introduce the spatial force correlation 
measurement which provides a natural definition for a long-range correlation length and in 
Section |V] we show how the size of the correlation length affects the contact forces between 
grains. 



II. GRANULAR SHEAR FLOW BASIC CONSIDERATIONS 

If no external force is applied to a dry granular material in the absence of gravity, it 
quickly loses all its kinetic energy in dissipative collisions, and each grain comes to rest. 
If this occurs for a dilute system there are no residual contacts between any grains and 
the total energy is zero. However, for granular materials with larger densities, there are 
contacts between grains in the relaxed state and a non-zero residual energy remains due to 
grain deformation and friction. 

If a shear stress is then applied to the system, motion only occurs if the stress is large 
enough to overcome the energy stored in the contacts. The minimum stress needed to initiate 
motion is called the yield stress, which is zero below a critical packing fraction v c and is an 
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increasing function of packing fraction above v c |35|, |36j, 1371 ] . 

For granular shear flows with v > u c , previous research has demonstrated that the stiffness 

n 

of the grains plays an important role at all values of the shear rate [38] . This is because grains 
are not able to rearrange to a configuration where no contacts exist and the system moves 



between different configurations where the grains are always de 
v > v c are characterized by slowly moving quasi-static flows 



35, 
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ormed. Shear flows with 



40| . where force balance 



371 ]. where motion ceases for stresses below the 



is upheld at all times, and jamming 
the yield stress. 

Conversely, for granular shear flows with v < u c , it has been demonstrated that the 
the grains can always be taken large enough so that it plays no role in the dy- 



stiffness o 
namics 38 
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43] . This is because grains are always able to rearrange to find free volume 



and the system moves between different configurations with minute grain deformation. In 
this regime inertial terms are dominant and an invariance in Newton's equations 44| shows 
that the dynamics are controlled exclusively by the shear rate 7 45[]. For these inertial flows 
the stress tensor is proportional to 7 2 , which is referred to as Bagnold's scaling and has been 
observed in experiments 46] and simulations 
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45]. 



In this paper we consider only inertial flows of granular materials with v < u c , in the 
limit of infinite stiffness. This allows us to investigate a wide range of packing fraction, from 
zero to u c , and explore how correlations between grains naturally emerge in this range. In 
particular, we measure spatial force correlations and observe that they decay exponentially 
with a characteristic length-scale £. By measuring distributions of contact forces, we demon- 
strate that only binary collisions are relevant for small £ and simultaneous contacts have 
quantitative effects on contact forces for large £. This allows us to define an important pack- 
ing fraction vi> c , where the crossover occurs between binary collisions and force networks. 
This crossover is accompanied by a signature in the contact force probability distribution 
function that has been observed previously in experiments on hopper flow [47 ] . 

The transition at u^ c is purely dynamical and would not occur if 7 = 0. We find that the 
value of Ub c depends on the amount of energy dissipated at each contact, but it is always the 
case that Z4, c is strictly less than v c . This separates the inertial regime into a dilute inertial 
regime where only binary collisions are relevant and a dense inertial regime where clusters 
of contacting grains form into complex force networks. The network transition at v^ c arises 
due to steric exclusion and is independent of the roughness of the grains. A schematic of 
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FIG. 1: Schematic of the regimes of granular shear flow. Here we explore inertial flows with v < v c 
and show that there is a network transition at v\, c where force chain networks begin to form and 
affect the dynamics. 

the different regimes available to a granular shear flow is displayed in Figure [TJ 
III. SIMULATIONS 



We perform simulations of inertial granular shear flow in two- dimensions using the Con- 
tact Dynamics (CD) algorithm. The CD algorithm was developed to model the dynamics 



of dense collections of rigid grains [48 
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51 
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53| and can include the effects of long 



lasting, simultaneous contacts between groups of grains. It is a fully dynamical and ather- 
mal algorithm, with the motion of each grain determined at every time step by integrating 
Newton's equations. 

To carry out the integration, the algorithm provides values for the forces between each 
pair of contacting grains. In dry granular materials, contact forces arise due to deformation 
of the grains upon contact and friction between grains. The CD algorithm considers the 
rigid limit in a self-consistent way where forces take on the precise value necessary to p revent 



deformation, uphold Coulomb friction, and adhere to Newton's equations [48j, |49j, [50|, [5l| . If 
we use & lJ to denote the unit vector connecting the centers of two contacting grains labeled i 
and j, then the deformation produces a normal force in the direction of a %3 and the interplay 
of friction and deformation produces a tangential force perpendicular to a 13 . These contact 
forces are dissipative and depend sensitively on the velocities of the colliding grains 54]. 
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The CD algorithm determines the forces arising from grain deformation by assuming 
that grains are infinitely rigid and setting constraints on the total energy dissipated in each 
contact. For two grains in contact with a relative velocity v /y , the algorithm determines a 
contact force such that the relative velocity in the next time step v*- 7 is given by 

v y . & ij = _ ev «j . & ij . y « x & ij = e yv x ^ 

In this way, the relative velocities are altered by restitution coefficients in the normal direc- 
tion (e) and tangential direction (et). This procedure ensures that the momentum transfered 
between grains with non-zero relative velocities is at least that which is supposed to arise 
from collision. Friction is included by assuming that the grains have a coefficient of friction 
/i. If the ratio of the tangential to the normal force exceeds /i, then the grains are allowed 
to slip with a tangential force equal to \x times the normal force. Using these dynamical 
constraints, contact forces can be determined at each time step. 

It is important to mention one subtlety of the algorithm: because a constant time step is 
utilized, many contacts can (and do) occur in each time step. In this case, all contacts are 
assumed to occur simultaneously and the dynamical rule in Equation ([T|), along with the fric- 
tion constraint, is applied to each contact. This leads to a set of coupled algebraic equations 
for each group of contacting grains that is solved using an iterative method. Therefore, the 
value of a specific contact force depends not only on properties of the contacting grains, but 
also on the properties of other grains in the connected cluster. Physically, this corresponds 
to the effects of contact forces propagating through a network of grains much faster than 
contacts are created and/or destroyed. The CD algorithm explores the limit of infinitely 
rigid grains by assuming that the propagation occurs instantaneously. 



The limit of infinitely rigid grains has also been explored in other simulations 



38 
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42 



431 ] . In one study 38j, Campbell simulated soft spheres with a rigidity k in shear flow at 
constant volume. He found that, for any packing fraction below a critical value u c , k can be 
taken large enough so that no macroscopic variables depend on its value. This study, along 



with others 
flow. 



41 



42j, |43J, establishes the validity of the "hard-sphere" limit in granular shear 



Shear flows of hard-sphere granular materials always occur in the inertial regime 44|, |45 |. 
Because the stiffness of the grains does not play a role, the interactions between grains are 
mediated by dimensionless quantities, as in Equation ([1]), and the only timescale is set by the 
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FIG. 2: Snapshot of a granular simple shear flow. Each grain has an average velocity in the 
x-direction given by jy, where 7 is the shear rate. The center of the cell is defined as x = y = 0. 



value of the external shear rate 7. Therefore, by dimensional analysis, it must be the case 
that the pressure and the shear stress are proportional to A f 2 . This is Bagnold's scaling 46], 
which is a central characteristic of granular flows in the inertial regime. Since hard-sphere 
granular flows must always occur in the inertial regime, our simulations are constrained 
to lie below the jamming transition at v c and always have zero yield stress. However, we 
can explore granular flows arbitrarily close to v c and we have conducted simulations up to 
packing fraction of 0.84. 

The value of v c in our system is ultimately determined by the distribution of grain sizes. 
In the simulations presented here we use a polydisperse collection of circular grains in two 
dimensions. The grain diameters are chosen randomly from a flat distribution with minimum 
and maximum diameters given by o ± A, with a = 1.4 and A = 0.26a. We have found that 
this amount of polydispersity restricts crystallization. We estimate that the value of u c for 
this grain polydispersity, given the shearing algorithm, is between 0.84 and 0.85. 

Finally we comment on the shearing algorithm. We use a numerical procedure to create a 
simple shear flow without any boundaries. This is done by using the Lees-Edwards boundary 
conditions j^j] along with the shod equations of motion [56] to initiate shear flow immedi- 



ately. We discuss these procedures at length in an earlier paper [45j. For our purposes here, 
it suffices to state that this procedure produces a translationally invariant simple shear flow 
with a linear velocity profile. There are no other gradients in the system and the stress 
tensor, granular temperature, and packing fraction are all constant. Although this is an 
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idealized simulation, we have shown that the resulting rheo 
found in steady state incline flow, far from the boundaries 
shear flow simulation is shown in Figure [2J 

IV. FORCE CORRELATIONS 

In this section we present the central measurement of the paper. We investigate spatial 
force-force correlations in steady state shear flows through the measurement of the correla- 
tion function 

C(£) = ^ =lZ ^ =1 „ -. (2) 

EiiF l -F* 

In this equation, F* is the total vector force (sum of contact forces) experienced by a grain 
% at position r* and the sums are taken over all grains that have at least one contact. The 
distance £ ranges over the entire system size and is not limited to grains in direct contact. 
We take the average value of C{£) over at least 5000 time steps in steady state shear flow. 
A non-zero value of C(£) reveals that, on the average, two grains separated by a distance £ 
have forces that are correlated. 

For perfectly rigid granular materials, where there is only a repulsive interaction between 
grains at contact, C (£) gives a quantitative measurement of the average effect of force chains 
of length £ in the material. A non-zero value of the correlation indicates that two grains 
a distance £ are connected through a cluster of simultaneously contacting grains and the 
force from one grain is being transmitted through the network to the other grain. It thereby 
establishes that simultaneous contacts exist and that forces propagate through networks. 
Positive values of the correlation correspond to situations where the total forces on each 
grain tend to be aligned. 

Because we make the measurement of C(£) while the material is in steady state shear 
flow, the correlations do not reveal the presence of static structure. Instead, because contacts 
between grains are being created and destroyed by the overall flow, the correlation function 
gives information on the average size of dynamic structures that are fluctuating in both 
space and time. 

We will demonstrate that the correlation function depends on the vector distance £ = i£ 
between pairs of grains, and that it decays exponentially with £. In the following sections, we 
first investigate the dependence of C (£) on the magnitude of £, thereby defining an isotropic 

8 



ogy is identical to the rheology 



45 1 . A screenshot of our simple 



correlation length £. Then we investigate the full dependence of C(£) and obtain the full 
functional form of £(#), which depends on the angle 9 between pairs of grains. 

A. Measurements of C(£) and the length scale £ 

We begin by measuring the isotropic part of the correlation C(£) in two dimensional 
simulations by averaging C(£) over all directions £. In Figure [3] we plot this directionally 
averaged correlation function C(£) for a frictionless material with e = 0.25. 
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FIG. 3: (a) The force correlation function C{€) for e = 0.25 and v = 0.81, where i is the distance 
between grains and a is the average grain diameter, (b) The logarithm of the magnitude, log |C(^)|, 
for packing fractions of v = 0.6, 0.77, and 0.81, illustrating the exponential decay of the correlations. 
The lines correspond to the function e - ^, where £ is determined from Equation (J3j) and plotted 
in Figured! 

The logarithm of the magnitude of the correlation function log | C (€) | is also plotted in 
Figure [3) From these plots we observe that the magnitude decreases exponentially, although 
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the decay is complicated by an oscillating function that accounts for the sign of C{£). The 
form of this oscillating function is not universal and depends on the exact value of density 
and restitution. Nevertheless, as a first approximation, we express the correlation function 
as C{£) ~ exp(— £/£), which introduces a length scale £ that quantifies the large I behavior 
of the correlations. The data in Figure [3] illustrates that this length scale increases with 
packing fraction. 

We find that the value of £ is well approximated by the equation 

? j™d£C(£) ' 1 J 

In Figure [3] we plot exp(— where £ is determined for each density from Equation ([3]), 
and we observe excellent agreement with the measured exponential decays of C (£) . Although 
Equation ([3]) has the same form as an expectation value in statistical physics, that is not 
the proper interpretation here. Rather, Equation ([3]) is simply a combination of integrals 
that gives the coefficient of any exponential function. 

In Figure H] we plot measurements of £, determined from Equation ([3]), for all of the 
frictionless granular flows we have simulated. The value of £ quantifies the average extent 
of force chains in the system and is an increasing function of density for each value of 
restitution. Higher values of restitution have lower values of £. This is consistent with 
dynamic force networks forming spontaneously as the material is sheared due to energy 
dissipation upon contact. Indeed, in the limit of e — > 1, the data in Figure H] suggests that 
no correlations exist. The formation of correlations is therefore related to the observation of 
inelastic collapse in sheared granular materials [5?]]. Collapsing grains contact increasingly 
often, allowing a single grain to have multiple contacts, even in the limit of perfectly rigid 



grains 



58]. 



For large values of packing fraction, the data for £ increases rapidly. The maximum 
packing fraction we are able to simulate is z/ max = 0.84, but we were not able to determine 
the value of £ at z/ max for all values of restitution. This is because the length scale of the 
correlations becomes larger than the maximum system size of 100 x 100 grains that we 
are able to simulate. For the high values of packing fraction, we conduct simulations for 
packing fractions in increments of 0.01. When e < 0.5 we reach a packing fraction where, if 
we increase the packing fraction by 0.01, £ is too large to be measured. The rapid growth 
of the correlations at high packing fraction suggests that £ is diverging. 
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FIG. 4: Main Figure: The length scale £ for frictionless granular shear flows, normalized by 
the average radius a, plotted for a wide range of packing fraction and restitution e. The value of 
£ increases rapidly for large packing fraction and asymptotes to £ei/°" — 0.785 for small packing 
fraction. Inset: Data for e = at multiple values of the time step. Here dt corresponds to the 
time step we use in this paper, which is smaller for larger packing fractions. The value of £ does 
not depend on the time step. 

For small values of packing fraction, Figure H] shows the value of £ approaching a limiting 
value of 0.785o\ This corresponds to the smallest correlation possible, which arises when 
the only relevant interactions are binary collisions between grains. For these small packing 
fractions, the sign of C (£) is exclusively negative and the exponential decay is not observed. 
Indeed, in the limited simulations we have conducted with a mono-disperse collection at 
very low packing fraction, we observe that C(£) is zero for all values of I except i = a. In 
the case of a polydisperse collection of grains, the form of C(£) for dilute flows depends on 
the relative probabilities to have binary interactions between grains of different sizes and 
exponential decay is not observed since the correlation function equals zero for all values of 
I larger than the maximum grain diameter. Even though exponential decay is not observed 
in the limiting case of binary collisions between grains, Equation ([3]) still provides a useful 
measure of correlation length that matches the exponential decay at higher densities. 

The limiting value of 0.785 is related to the probability to have a binary collision between 
grains of different diameters and is equal to the correlation length that would be expected in 
a perfectly elastic system where e = 1 and no energy is dissipated upon interactions between 
grains. Thus we denote the limiting value as £ e i and expect that its numerical value is related 
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to the distribution of grain sizes. We have conducted a limited number of simulations to 
determine the behavior of £ e i for different grain distributions, and the results are plotted 
in Figure Each grain distribution we consider is a flat distribution with minimum and 
maximum grain diameter given by a ± A. As expected, £ e i/cr = 1 for A = 0. For large 
values of A, £ e i depends linearly on A since the largest grains set the correlations. The grain 
distribution we consider in this paper corresponds to A = 0.26cr and is designated by a solid 
data point. Although this grain distribution happens to occur close to the minimum, this 
is not relevant to the arguments in this paper. We present the data in Figure [5] to better 
understand the origin of £ e i, and in particular why it is not always equal to one. We will see 
that the ratio £/£ e i seems to signal the transition between different regimes of shear flow, 
and it asymptotes to one for small packing by definition. 




0.1 0.2 0.3 0.4 0.5 

A/a 



FIG. 5: The equilibrium correlation length £ e i as a function of the grain distribution of the system. 
We consider flat distributions with minimum and maximum grain diameter given by a ± A. The 
distribution considered in this paper is indicated by a solid circle. 

The data we have presented for frictionless grains can also be extended to situations where 
fi > 0. Adding friction introduces a non-zero tangential force at each contact that tends 
to increase geometrical frustration. Measurements of £ for frictional systems are plotted in 
Figure [6] for e = and three friction coefficients fi. We see that the value of £ diverges sooner 
for systems with friction. This is to be expected, since the additional constraints associated 
with the tangential forces can only increase the correlations. Interestingly, the value of £ for 
different friction coefficients is not different below a rather large packing fraction. This is 
because clusters of simultaneously contacting grains must exist before the sticking effects of 
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friction alter the dynamics and correlations. From a qualitative viewpoint, the behavior of 
£ for frictional systems is analogous to the frictionless case- the length-scale is observed to 
diverge at finite packing fraction and asymptote to £ e i at low packing. 
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FIG. 6: The length scale £, normalized by average diameter a, for granular shear flows with friction 
coefficient /i and e = 0. Larger values of /i correspond to larger £. The value of £/<r is always larger 
than £ el /cr = 0.785 




B. Connection to Jamming 



The data in Figures 0] and [6] suggests that £ diverges at a finite packing fraction u c that 
may depend on the friction coefficient and restitution coefficient. This divergence is related 
to the jamming transition 



59j ] in granular materials, where the shear modulus becomes 
non-zero and the system is able to sustain a shear stress without yielding. In order to 
make the transition from a flowing shear state to a jammed state, it is necessary that a 
correlation length approach the size of the system, and diverge in the thermodynamic limit. 
This is because force chains must percolate from the upper to lower shearing wall in order 
to counteract the applied shearing force. The correlation length £ quantifies the notion of 
force chains and we expect that the observed divergence of £ is a necessary condition for 
jamming. 

However, there is no guarantee that the divergence of £ is also a sufficient condition 



for jamming, since it is possible that force chains percolate 



Nevertheless, both theories 
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ong before the system jams. 
61| and simulations [35|, |36j, |37J] have found that percolation 
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and jamming occur simultaneously, which suggests that a granular system jams if, and only 
if, £ diverges. If this holds, then we would expect the jamming transition to occur at lower 
packing fraction as the friction between grains increases. 

C. Anisotropy in the Angular Dependence of £ 

We have also conducted measurements of the full directional dependence of C(£). In this 
case, we observe that the decay of the correlations can still be described by an exponential, 
but the value of the correlation length depends on the orientation £. In two dimensions this 
orientation can be quantified by the angle 9 between £ and the x-axis. In Figure [7J we plot 
the angular dependence of £ for frictionless flows with e = and three different packing 
fractions. These plots reveal that £ is not isotropic and the correlation between a pair of 
grains depends on their orientation. For a simple shear flow, correlations are created along 
the compressional axis due to a higher number of grain interactions, and is destroyed along 
the dilational axis due to contacts being lost. 




FIG. 7: The angular dependence of the length scale £(0) for a frictionless shear flow with e = 
and three different packing fractions v. The length scale is greatest along the compressional axis 
of the shear flow. 

We note that the maximum value of the correlation length occurs at approximately the 
same angle for each packing fraction in Figure [7J This trend is followed for other packing 
fractions and restitution coefficients as well. In Figure [HI we plot the angular dependence of 
the length scale divided by its average value: £(#)/£. We notice that this collapses the data 
for a large range of packing fractions and restitution coefficients onto one curve. The collapse 
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is not perfect, especially along the dilational axis of the shear flow where the correlations 
are small. We believe this to be due to the small number of collisions that occur on this 
axis, which makes gathering statistics difficult. 
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FIG. 8: The normalized angular dependence of the length scale, for frictionless shear flows 

with different restitution coefficients and packing fractions. The data is well characterized by 
Equation (jU), plotted as a line. 



The common collapsed curve for all of the data in Figure M shows that £ is anisotropic. 
The angular dependence of £ can be approximated as a Fourier series, only including terms 



that are 7r-periodic. We find that, as is the case for other anisotropies in granular flow 
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64 . l65l . l66l . l67j ] , the functional form of £ (9) is well characterized by the first two terms in 



the Fourier series 

£(0) = £(l-a o sm2[0-0o]), ( 4 ) 
where £ is the average value of the correlation length. The solid curve in Figure [8] is a 
fit to this equation, which estimates the parameter values as = 0.21 and 9q = 0.013. 
The value of 9q is consistent with the axis of maximum compression for the data we have 
gathered. This implies that larger values of £ occur near 9 = 3tt/4 because the compression 
induced by the shear flow causes more grains to come into contact. Near 9 = it/ 4, where £ 
is a minimum, dilation reduces the magnitude of correlations and thereby the length scale. 
However, we have found no simple explanation for the value of a®. We suspect that the 



numerical value of a® should be related to the anisotropy in the contact distribution [65 ]. 
However, because the anisotropy in £ remains constant over a range of packing fractions and 
restitution coefficients where the anisotropy in the contact distribution is not constant, the 
value of a may have a simpler origin. 
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V. THE EFFECTS OF CORRELATIONS ON CONTACT FORCES 



In the previous section we presented measurements of a correlation length £ that diverges 
at the jamming transition and asymptotes to an elastic value £ e j at small packing fraction. 
This length scale captures the decay of force correlations and is related to the emergence of 
clusters of simultaneously contacting grains. In very dilute systems only binary collisions 
are relevant, £ = £ e i, and contact forces can be determined from the properties of the two 
colliding grains. However, as the packing fraction is increased, £ also increases and the 
contact force between any two grains will depend on properties of the other grains in the 
cluster. This is because, in the rigid grain limit, forces propagate instantaneously through 
the network. Because the growth of £ is closely related to the nature of the force transfer- 
binary collisions or force networks- we expect the contact forces to depend on the value 
of £. A useful way to explore properties of contact forces is to measure the contact force 
probability distribution function (pdf) P(F). This function encodes the statistics of the 
contact forces: P(F)dF is proportional to the number of contact forces in the range F to 



To make the connection between contact forces and long range spatial force correlations 
explicit, we begin by demonstrating that the contact forces between pairs of grains can not 
be determined simply by assuming binary collisions when £ is large. To illustrate this point, 
we compare the statistics of the actual contact forces to the forces we would calculate if we 
assumed that only binary collisions occurred. If we make the binary collision assumption, 
then the dynamical rule in Equation (pQ), along with momentum conservation, determines 
the normal impulse in each collision. Dividing this impulse by the algorithm time step yields 
the average force that would arise over the time step dt. We label this force F^ J C , where i and 
j represent the colliding grains and the label "be" reminds us that this force only applies to 
purely binary collisions. It is simple to show that the value of the binary force is given by 



where e is the normal restitution coefficient, \x = m l m? j[m % + m 3 ) is the reduced mass, v n 
the pre-collisional velocity of grain i, and a 13 is the unit vector connecting the centers of 
grains i and j. All of these terms can be measured directly from simulations. 

In Figure [9] we plot measurements of the contact force distribution function P(F) for six 
different values of restitution and packing fraction. In each figure we compare P(F) with 



F + dF. 




(5) 
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FIG. 9: The data points are P(f) for systems with e = 0.75 and growing values of £/£ e i> where 
/ is the contact force F divided by the average value (F). This data is compared with the line, 
where the force is determined from Equation ([5]) which assumes that only binary collisions occur. 
We have normalized both the total forces and binary collision forces in each plot by their average 
values. There is excellent agreement for £/£ c i < 1-25. For larger values of the correlation length, 
clusters of interacting grains form and the binary collision assumption does not fit the data. 
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the statistics of the binary forces P(Fb c ). If these two functions are the same then contact 
forces are well approximated by only considering binary collisions; if the functions differ 
then we know that clusters of contacting grains affect contact forces. 

We indicate the value of £/ £ e i for each plot in Figure M and immediately see that for 
small values of £/£ e i the data for P(F) is well fit by the line, which is a measurement 
of P(-Fb c ). However, as £ increases, the presence of force networks changes the nature of 
the contact forces and we can no longer make accurate predictions by assuming that only 
binary collisions occur. The value £/£ e i ~ 1-25 serves as a rough upper bound for the 
regime where the binary collision assumption is reasonable. Similar values of £/£ e i signal 
the transition between collisional and non-collisional flows for all restitution coefficients and 
packing fractions we have investigated. 

This is not surprising since £/£ e i > 1-25 comprises a region where force networks have 
formed and simultaneous contacts occur. In this regime, in order to calculate the force 
between two grains, it is not sufficient to only consider the properties of the two contacting 
grains. Rather, all of the grains connected in the force network play an important role. This 
is because the two contacting grains are being pushed together by the other grains in the 
cluster and the contact force is equal to the binary collision contribution from Equation (jSJ) 
plus a contribution from the cluster. 

We conclude from Figure [9] that £/£ e i ~ 1-25 separates the region where only binary 
collisions occur from the region where force networks begin to form and affect contact forces. 
This is an approximate value of the length scale at which there is a clear deviation between 
binary collisional forces and total forces. Measurements presented later in this paper, and in 
the companion paper [^(J, also place the network transition at £/£ e i ~ 1-25. Although this 
value of the length scale does not change with restitution coefficient or friction, we would 
expect it to depend on the grain size distribution. 

The techniques we used to determine the crossover in Figure [9] can only be used in sim- 
ulations where the position, velocity, and force on every grain is always known. However, 
we have also found a signature of the transition that can (and has been) observed in exper- 
iments of granular flows. This signature relates to the small force behavior of the contact 
force distribution function. 

In Figure [10] we present data for the contact force distribution function P(F). In par- 
ticular, we plot log P(f) for many different values of the restitution coefficient and packing 
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FIG. 10: Measurements of the contact force probability distribution function P(f), where / is the 
value of the contact force divided by the average contact force. The different curves have been 
offset so each curve is easy to see. (a) Data for v = 0.7 and different values of the restitution 
coefficient e. (b) Data for e = 0.75 and increasing values of packing fraction v. The curves in each 
plot are also labeled by their associated value of £/£ e i, and we observe that the peak is present in 
P(f) only if e/Cel < 1-28. 

fraction, where / is equal to the normal contact force F divided by the average normal force 
(F) . All of these curves correspond to frictionless materials, but we have observed that the 
statistics of the normal forces display the same behavior for frictional systems. Our mea- 
surements of P(f) have been averaged over 5000 time steps in steady state shear flow, and 
the different curves are vertically displaced in the figure so each can be clearly seen. Each 
curve is labeled by the value of £/£ e i and we immediately see that the behavior at small / 
depends on the value of £. For £/£ e i < 1-25 there is a clear peak, whereas for £/£ e i > 1.25, 
the peak disappears and the maximum occurs at / = 0. 

This measurement once again defines a crossover around £/£ cl = 1.25. This is the tran- 
sition where the microscopic interactions change from being dominated by binary collisions 
between grains to being dominated by clusters of grains, forming force networks of size £/£ e i- 
When this network transition occurs, the peak disappears and the most likely force is no 
longer equal to the average force. This is because grains have spontaneously formed into 
transient clusters and the greatest number of contacts are simply rolling over each other, 
which produces a very small normal force. This moves the peak to / = once the transi- 
tion has fully developed and the average force is not representative of most of the forces. 
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Additionally, as force-networks become long-ranged, the data shows that there is a greater 
probability of large forces, which arise from a large number of grains in a cluster compressing 
two contacting grains. 

We have also taken data of P(F) close to the network transition point for granular flows 
with e = 0.75. In Figure [TT] we plot the value of the force -F max for which P(F max ) is 
a maximum, as a function of packing fraction. Between packing fractions of 0.766 and 
0.767, our data shows that the value of / max (equal to F max divided by the average force 
(F)) makes a jump from approximately 0.2 to zero. Within this small range of packing 
fraction, the contact force pdf loses its peak. The abruptness of the transition suggests 
that different physical processes are occurring on either side of the transition and that the 
network transition may be quite sharp. 
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FIG. 11: The value /max at which the contact force pdf P(f max ) is a maximum, as a function of 
packing fraction, for e = 0.75. The maximum of P(f) occurs at positive / for v < 0.766 and at 
/ = for v > 0.767. 



The signature of the transition evident in our measurements of P(f) has also been ob- 
served in other simulations and experiments, but has never been connected to the formation 



of large scale structure. Simulations of granular hopper flow |68|, |69| , conducted using Event 
Driven simulations where only binary collisions are allowed to occur, have reported that as 
the hopper aperture is reduced and the density of the packing increases, P(f) begins to lose 
its peak. This is consistent with our results from Figure [10] and suggests that the correla- 
tion length £ is relevant in more than just shear flows of granular media. Additionally, the 
same behavior in P(f) has been observed in experiments on hopper flow 47|, which lends 



20 



credibility to the result and suggests that it is not an artifact of the simulation methods 



used here and in Refs. 



68 



691 ]. but rather a real effect in granular flows. 



The results we have cited from hopper flow, and others 70j, have been used to challenge 
the belief that the formation of a peak in P(f) is a signature of the jamming transition 71, 
Q. In a wide variety of contexts, including incline flow |t1 | . quasi-static flow 22, 72], 
and jammed granular materials 3l|, |32|, |33|, |73| it has been observed that P(f) exhibits a 



maximum at / = (no peak) if the system is flowing, while a peak at non-zero / forms 
as the systems jams. These observations, coupled with similar results in Lennard- Jones 
glasses and foams [iJ], have been used to bolster the claim that the formation of a peak in 
P(f) is a generic characteristic of the jamming transition, and a necessary condition for the 
appearance of a yield stress. 

Our observations reveal that, in fact, there are two important transitions encoded in the 
small / behavior of P{f). First, at a low packing fraction, there is the network transition 
where interactions between grains change from binary collisions to force networks. This 
occurs in the inertial flow regime and is accompanied by a change in P(f) where the peak 
that was present for small densities disappears and the maximum value of P(f) occurs at 
/ = 0. Then, as shown elsewhere, there is another transition at higher packing fraction 
where the system develops a yield stress and the peak reappears in P(f). 

In summary, we have measured contact force distribution in this section to determine 
the effects of long-range correlations. We find that £/£ e i = 1-25 separates the regime where 
only binary collisions occur from the regime where force networks form. This observation 
allows us to split the inertial regime where hard-sphere granular flows exist into two distinct 
regions. At low packing fraction there is a "dilute regime" where binary collisions are the 
dominant microscopic interaction and £/£ e i < 1-25. At high packing fraction there is a 
"dense regime" where force networks exist but do not percolate through the system. This 
dense regime is characterized by clusters of interacting grains with an average extent £ / £ e i in 
the range of 1.25 < £/£ e i < oo. For dilute flows where only binary collisions occur, a peak is 
visible in P(f); as force networks begin to appear in the dense regime, the peak disappears. 

The crossover at £/£ e i = 1-25 defines the transition between interactions dominated by 
binary collisions and interactions dominated by force networks. Therefore, we use this value 
to define Z4, c (e), which is the value of the packing fraction below which only binary collisions 
are relevant (shown schematically in Figured]). This function is plotted in Figure fT2l using 
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FIG. 12: The value of v\> c as a function of restitution coefficient e. Below i^ c (e) only binary 
collisions are relevant and above fb c force networks emerge. As e approaches unity the network 
transition is concurrent with the jamming transition at v c ~ 0.845. 

the data for £ from Figure HI This plot is for frictionless materials fi = 0, but increasing the 
value of /i does not change the curve. This is because, as we saw in Figure El the effects 
of friction do not take hold until £ is much larger than 1.25£bc- Thus is an important 
packing fraction for all values of friction and we see that it is an increasing function of e. 
This is because larger e produces less energy dissipation and restricts grain clustering. The 
data in Figure [12] implies that as e — > 1 then v^ c — » v c , which is the packing fraction at 
which the system jams. This means that as grains become perfectly elastic and no energy 
is dissipated at contacts, then binary collisions describe the interactions for all values of 
packing fraction in the inertial regime. 

VI. CONCLUSION 

We have measured correlations between grain-forces in inertial shear flows. These cor- 
relations are long-ranged, decaying with a characteristic length scale £ that diverges at the 
jamming transition and asymptotes to an elastic value £ e i m the dilute limit. By investigat- 
ing the statistics of contact forces between grains, we have shown that £/£ e i — 1-25 splits the 
inertial regime into dilute flows, where all forces arise from binary collisions between grains, 
and dense flows, where force chain networks begin to form. We denoted as the packing 
fraction at which £/£ e i = 1-25 and found that the value of v\> c depends on the restitution 
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FIG. 13: State diagram of granular shear flow, plotted as a function of the packing fraction v and 
the restitution coefficient e. For our system, we estimate that u c ~ 0.845, which is represented by 
the vertical line on the plot. The value of v c separates the inertial regime from the quasi-static 
regime. In this paper we have focused on the inertial regime and shown that there is an important 
packing fraction, z^, c , where £/£ e i = 1-25. Our measured values of u^ c depend on the restitution 
coefficient and are plotted on the diagram. For v < the flow is in the Dilute Inertial Regime 
where microscopic interactions consist of binary collisions; for v > the flow is in the Dense 
Inertial Regime where long-ranged force networks begin to form. 



coefficient, but is always less than the packing fraction u c at which the system jams. This 
phenomenology is illustrated in Figure [T3J 

The crossover from the dilute to dense regime is accompanied by a qualitative change in 
the nature of contact forces between grains, measured using the contact force distribution 
function P(f). For shear flows in the dilute regime P(f) has a peak at non-zero /, whereas 
for shear flows in the dense regime there is no peak. This raises interesting questions about 
the nature of the transition that occurs at u\, c . We observe that the average coordination 
number z is equal to one for flows with v < v\, c and is greater than one for v > v\, c . Further 
investigation of the behavior of granular flows very close to is needed to bring further 
insight to the physics near the transition. 

From a wider viewpoint, the presence of long range correlations changes the assumptions 
that can be made when modeling properties of the system. We explore the theoretical 
implications of correlations in the companion paper 20] . 
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